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Abstract 

In the mathematical modelling of sediment com- 
paction and porous media how, the rheological 
behaviour of sediments is typically modelled in 
terms of a nonlinear relationship between effective 
pressure p e and porosity <fi, that is p e = p e (<fi)- 
The compaction law is essentially a poroelastic 
one. However, viscous compaction due to pressure 
solution becomes important at larger depths and 
causes this relationship to become more akin to 
a viscous rheology. A generalised viscoelastic 
compaction model of Maxwell type is formulated, 
and different styles of nonlinear behaviour are 
asymptotically analysed and compared in this 
paper. 

Citation detail: X. S. Yang, Nonlinear viscoelastic 
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1 Introduction 

Drilling mud is used in well-bores drilled for oil ex- 
ploration to maintain the integrity and safety of the 
hole. The mud density must be sufficient to prevent 
collapse of the hole, but not so high that hydrofrac- 
turing of the surrounding rock occurs. Both these 
effects depend on the pore fluid pressure in the rock, 
and drilling problems occur in regions where abnor- 
mal pore pressure or overpressuring occurs, such as 
in the sedimentary basins of the North Sea, where 
pore pressure increases downward faster than hy- 
drostatic pressure. Such overpressuring can affect 
oil-drilling rates substantially and even cause se- 
rious blowouts during drilling. Therefore, an in- 
dustrially important objective is to predict over- 
pressuring before drilling and to identify its pre- 
cursors during drilling. Another related objective 
is to predict reservoir quality and hydrocarbon mi- 
gration. An essential step in the achievement of 
such objectives is the scientific understanding of 
their mechanisms and the evolutionary history of 
post-depositional sediments such as shales. 



Shales and other fine-grained compressible rocks 
are considered to be the source rocks for much 
petroleum found in sandstones and carbonates. At 
deposition, sediments such as shales and sands typ- 
ically have porosities of about 0.5 ~ 0.75 (Lerche, 
1990). When sediments are drilled at a depth of say 
5000 m, porosities are typically 0.05 ~ 0.2. Thus 
an enormous amount of water has escaped from the 
sediments during their deposition and later evolu- 
tion. Because of the fluid escape, the grain-to-grain 
contact pressure must increase to support the over- 
lying sediment weight. Dynamical fluid escape de- 
pends lithologically on the permeability behavior of 
the evolving sediments. As fluid escape proceeds, 
the porosity decreases, so the permeability becomes 
smaller, leading to an ever-increasing delay in ex- 
tracting the residual fluids. The addition of over- 
burden sediments is compensated by an increase 
of excess pressure in the retained fluids. Thus 
overpressure develops from such a non- equilibrium 
compaction environment (Audet and Fowler, 1992; 
Fowler and Yang, 1998). A rapidly accumulat- 
ing basin is unable to expel pore fluids sufficiently 
rapidly due to the weight of overburden rock. The 
development of overpressuring retards compaction, 
resulting in a higher porosity, a higher permeabil- 
ity and a higher thermal conductivity than are nor- 
mal for a given depth, which in turn changes the 
structural and stratigraphic shaping of sedimentary 
units and provides a potential for hydrocarbon mi- 
gration. Therefore, the determination of the mech- 
anism of the dynamical evolution of escape of fluids 
and the timing of oil and gas migration out of such 
fine-grained rocks is a major problem. The funda- 
mental understanding of mechanical and physico- 
chemical properties of these rocks in the earth's 
crust has important applications in petrology, sed- 
imentology, soil mechanics, oil and gas engineering 
and other geophysical research areas. 

Compaction is the process of volume reduction 
via pore-water expulsion within sediments due to 
the increasing weight of overburden load. The re- 
quirement of its occurrence is not only the applica- 
tion of an overburden load but also the expulsion of 
pore water. The extent of compaction is strongly 



influenced by the burial history and the lithology 
of sediments. The freshly deposited loosely packed 
sediments tend to evolve as an open system, toward 
a closely packed grain framework during the ini- 
tial stages of burial compaction and this is accom- 
plished by the processes of grain slippage, rotation, 
bending and brittle fracturing. Such reorientation 
processes are collectively referred to as mechanical 
compaction, which generally takes place in the first 
1-2 km of burial. After this initial porosity loss, 
further porosity reduction is accomplished by the 
process of chemical compaction such as pressure so- 
lution (Fowler and Yang, 1999). 

Despite the importance of compaction and di- 
agenesis for geological problems, the literature of 
quantitative modelling is not a huge one, though 
the processes have received much attention in the 
literature, and the mechanism leading to pressure 
solution is still poorly understood. The effect of 
gravitational compaction was reviewed by Hcd- 
berg (1936) who suggested that an interdisciplinary 
study involving soil mechanics, geochemistry, geo- 
physics and geology is needed for a full understand- 
ing of the gravitational compaction process. More 
comprehensive and recent reviews on the subject 
of compaction of argillaceous sediments were made 
by Rieke and Chilingarian (1974) and Fowler and 
Yang (1998). 

Compaction is a density-driven flow in a porous 
medium, a fascinating multidisciplinary topic that 
has attracted attention from scientists with differ- 
ent expertise for a long time. Holzbecher (1998) 
provides an up-to-date comprehensive review of 
the previous work and state-of-art numerical meth- 
ods and softwares for modeling density-driven flow 
and transport in porous media where the constant 
porosity is used. Here, we will mainly model how 
porosity changes with time and depth, rather than 
using a constant density; thus an appropriate com- 
paction relation is vitally important. 

Nonlinear compaction models have been formu- 
lated in two ways. The early and most models used 
elastic or poroelastic rheology, and the compaction 
relation is of Athy's type p e — Pe{4>) (Gibson, Eng- 
land & Hussey, 1967; Smith, 1971; Sharp, 1976; 
Wangen, 1992; Audet and Fowler, 1992; Fowler 
and Yang, 1998). More recent models use a viscous 
rheology with a compaction relation p e = — £V.u s 
(Angevine and Turcotte, 1983; Birchwood and Tur- 
cotte, 1994; Fowler and Yang, 1999). The poroe- 
lastic models are valid for mechanical compaction 
while the viscous models mainly describe the chem- 
ical compaction such as pressure solutions. More 
recently, efforts have been made to give a more re- 



alistic visco-elastic model (Revil, 1999; Fowler and 
Yang, 1999). Fowler and Yang (1999) use a viscous 
law p e = — £V.u s to model compaction due to pres- 
sure solution, while Revil (1999) uses a poro-visco- 
plastic model, with a relationship between poros- 
ity strain and effective stress, to study pressure so- 
lution mechanism and its applications. However, 
there is no viscoelastic model which has been for- 
mulated to analyse the compaction problem on a 
basin scale, and most of the conventional studies 
are mainly numerical simulations. Obviously more 
work has yet to be done. This paper aims at provid- 
ing a unified approach to the compaction relation 
by using a visco-poroelastic relation of Maxwell 
type. The nonlinear partial differential equations 
are then analysed by using asymptotic methods and 
the analytical solutions are compared with numer- 
ical simulations. 



2 Mathematical Model 

For the convenience of investigating the effect 
of sediment compaction, we will assume a single 
species only. The sediments act as a compress- 
ible porous matrix, so that mass conservation of 
pore fluid together with Darcy's law leads to the 
equations: 

Conservation of mass, 

|(l-0) + V-[(W)u s ]=O, (i) 

^+V.(0u') = O, (2) 

Darcy's law, 

4>{u l -u') = --{W + Pl93), (3) 
A 4 

Force balance, 

V<r e - V[p ; ] - P<?j = 0, (4) 

where u' and u s are the velocities of the fluid and 
solid matrix, k and [i are the matrix permeability 
and the liquid viscosity, pi and p s are the densi- 
ties of the fluid and solid matrix, <x e is the effec- 
tive stress, p e is the effective pressure, j is the unit 
vector pointing vertically upwards, p l is the pore 
pressure, and g is the gravitational acceleration. In 
addition, a rheological compaction law is needed to 
complete this model. 



2.1 Poroelasticity and Viscous Com- 
paction 

The compaction law is a relationship between ef- 
fective pressure p e and strain rate e or porosity <f>. 
The common approach in soil mechanics and sed- 
iment compaction is to model this generally non- 
linear behaviour as poroelastic, that is to say, a 
relationship of Athy's law of the form p e = p e (4>), 
which is derived by fitting real data for sediments. 
Athy's poroelasticity law is also a simplified form of 
Critical State Theory (Schofield and Wroth, 1986). 
A common relation describing poroelasticity is 



and equation (12) can be rewritten as 



(5) 
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combining with the previous equation, we have 

Pe =Pe{<t>), 



(7) 



which is Athy's law for poroelasticity. A typical 
form of this constitutive relation (Smith, 1971; Au- 
det and Fowler, 1992; Fowler and Yang, 1998) is 



p e = ln(0o/<£) - (0o - 4>)- 



(8) 



However, this poroelastic compaction law is only 
valid for sediment compaction in the upper and 
shallow region, where compaction occurs due to 
purely mechanical movements such as grain sliding 
and packing rearrangement. In the deeper region, 
mechanical compaction is gradually replaced by 
chemical compaction due to stress-enhanced flow 
along grain boundaries from the grain contact ar- 
eas to the free pore, where the pressure is essen- 
tially the pore pressure. A typical process of such 
chemical compaction in sediment is pressure solu- 
tion, whose rheological behavior is usually viscous, 
so that it is sometimes called viscous pressure so- 
lution or viscous creep. 

The mathematical formulation of compaction 
laws for pressure solution is to derive the creep rate 
in terms of concentrations, grain size and geometry 
(usually spherical or cylindrical packings) , effective 
stress, grain boundary diffusion. This allows us to 
include the detailed reaction-transport process in a 
relation between strain rate and effective stress, al- 
though further simplifications are usually assumed 
such as steady-state dissolution and local reprecip- 
itation along the grain boundary. Rutter's creep 



law (1976) is widely used 



A k c WDgb 
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(7, 



(9) 



where a is the effective normal stress across the 
grain contacts, is a constant, cq is the equilib- 
rium concentration (of quartz) in pore fluid, p, d 
are the density and (averaged) grain diameter (of 
quartz). D g i, is the diffusivity of the solute in water 
along grain boundaries with a thickness w. Note 
that p e = —a and ikk = V • u s . With this, (9) 
becomes the following compaction law 



Pe 



-£V.u s . 



(10) 



This was first used by Birchwood and Turcotte 
(1994) to study pressure solution in sedimentary 
basins by presenting a unified approach to geopres- 
suring, low permeability zone formation and sec- 
ondary porosity generation. 

2.2 1-D Viscoelastic Compaction 

Following the discussions of elastic compaction 
(Fowler and Yang 1998) and viscous compaction 
(Fowler and Yang, 1999), we can generalise the 
above relations to a viscoelastic compaction law of 
Maxwell type 



V.ir 5 = 



1 Dpe _ 1 

K Dt tf e 



(11) 



Equivalently, we would anticipate that a viscoelas- 
tic rheology holds for the medium, involving mate- 
rial derivatives of tensors, and some care is needed 
to ensure that the resulting stress-strain relation 
be invariant under the coordinate transformation. 
This is not always guaranteed due to the com- 
plexity of the rheological relations (Bird, Arm- 
strong & Hassagcr 1977). Fortunately, for one- 
dimcnsional flow, which is always irrotational, the 
equation is invariant and all the different equations 
in corotational and codeformational frames degen- 
erate into the same form (Yang, 1997). In the one- 
dimensional case we discuss below, we can take this 
for granted. The 1-D model equations become 



dt 



0, 



d(f> d(4>u l 



(f>{u l -u s ) = 



dt 



dz 

,dp e 



= 0, 



(12) 



(13) 



du s 1 Dpe 
~~dz ~ ~K~Dt 



-G^-(p s - Pl )(l-cP)g}, (14) 
1 D 9 J 

-\ v - m = dt + u d- z ' (15) 



where G = 1 + 4?y/3£ is a constant describing the 
properties of the sediments, and r\ is the viscosity of 
the medium. By assuming that the densities p s and 
pi are constants, we can see that only the density 
difference p s — Pi is important to the flow evolution. 
Thus, the compactional flow is essentially density- 
driven flow in a porous medium (Holzbecher, 1 998). 

3 Non-dimensionalization 

Let the length-scale d is defined by 



d = { 



£rh s G y /2 



[Ps- Pi)g' 



(16) 



so that the dimensionless pressure is p = Gp e /(p s — 
pi)gd = 0(1). Meanwhile, we scale z by d, u s by 
tyl s , t by d/m s , permeability k by k . By writing 
k{(f>) = fc fc*, z = dz* , and dropping the aster- 
isks, we have 



-£ + £i<i-**i-o, 

d<t> | d{<jm l ) _ ^ 
dt dz 

J>(u l -u s ) = \k( ( j>)l-^-(l-4 



(17) 

(18) 
(19) 



du s 



D P D d s d 
2 ^-P, — = s — , (20) 



<9z (f - 0) 2 " Dt dt dz 
where 

MPs - Pi) 9 



A 



(21) 



In the above derivation, we have used the require- 
ment that the poroelastic case (8) result in the limit 
as the viscous rheology vanishes. 

Adding (17) and (18) and integrating, we have 



u s = -<j){u l - u s ), 



where u — <f)(u l — u s ) is the Darcy flow velocity. By 
using the equation (6), we have 



dp 



(23) 



»' = «p m [-f z -{!-*)]■ (24) 



1 D(j) 



Dp 



(1 -(b) Dt (1 - 0) 2 Dt P ' (25) 

The constitutive relation for permeability fc(0) is 
nonlinear, and of the form, 



fc(0) = (-f) m , m = 8. (26) 

00 



The boundary conditions at z = are 

- (1 - 0) = (or equivalents, u s = 0), (27) 



to, P = 0, 



(28) 



dp 



h = m(t) + \(^) m [ 7 f-(l -0)] at z = h(t), 
0o Cz 

(29) 

which is a moving boundary problem. 

We estimate these parameters by using values 
taken from observations. From the typical values 



of pi ~ 10 3 kgm~ 3 , p s ~ 2.5 x i0 3 kgm~ 



fc 



10~ 15 — 10- 20 m 2 , M - 10- 3 Nsm 2 , £ - 1 x 10 21 N 
s m~ 2 , m s ~ 300m Ma -1 = 1 x 10 -11 m s _1 , g ~ 
10ms- 2 , G ~ 1, d ~ 1000m we get A w 0.01 ~ 
1000. The only parameter A, which governs the 
evolution of the fluid flow and porosity in sedimen- 
tary basins, is the ratio between the permeability 
and the sedimentation rate. 



4 Asymptotic Analysis 

Since the nondimensional parameter A m 0.01 ~ 
1000 that controls the compaction process varies 
greatly, we can expect that the two limits A« 1 
and A> 1 will have very different features of poros- 
ity and flow evolution. In fact, A = 1 defines a tran- 
sition between slow compaction (A << 1) and fast 
compaction (A >> 1). We follow the asymptotic 
analysis of Fowler and Yang (1998) and Fowler and 
Yang (1999) to obtain some analytical asymptotic 
solutions. 

4.1 Slow Compaction (A <C 1) 



IfA<Cl,2~l,i~l,y~l, then u s < 1 and 

dtp 



(22) d± 



-gt ~ 0, and it follows that w O and D<j)/Dt 
f . Thus 



— « -A(1-0 O )^2, 



dt 



dp 



u s ^x[-f z -(i-4> )}, 

1 dej) 0o <9p 



(1 - O ) at (1 - 0o) 2 at 

which can be rewritten approximately as 
dp N/ 9 2 p , (1-0 O ) 2 w (l-0o) 2 A 



(30) 
(31) 

-P, (32) 



9t <9z 2 0o 



-p, A' = 



with the boundary conditions 

— w 1 - 0o, on z = 0, 



(33) 



(34) 



L 
0.2 




Figure 1: Comparison of numerical solutions (solid 
curves) with asymptotic solution (39) for A = 0.01 
at t = 3, 5. Where Z = z/h(t) is the scaled height. 



oo, 



(35) 



This is in fact equivalent to the case of conduction 
in a semi-infinite space with a constant flux at z = 
0. The solution in this case can be approximately 
expressed as 

(36) 
(37) 



(l-0o)V4A^ ierfc(C) + v / A 7 0o exp[ 



where 



/4A't 



and 



icrfc(C) = -^e"^ - Ccrfc(C). (38) 
This gives the approximate solution of as 



(i - <M te 



Ri O - o v4Vi ierfc(C) 

I j - O ] 

(39) 

thus compaction essentially occurs in a boundary 
layer near the bottom with a thickness of the order 
of a/A 7 - The comparison of this approximate solu- 
tion (39) (dashed curves) with numerical solutions 
(solid curves) is shown in Figure 1 for the values of 
A = 0.01,t = 5. The approximate solution is accu- 
rate when (0 j '(/)o) m C 1 so that t ~ \jm\f\ ~ 5, 
due to the fact that <p(z = 0) = 1 - 0(V\t). The 
agreement is clearly shown in the figure. 



4.2 Fast Compaction (A > 1) 

Either viscous or poroclastic fast compaction, is 
more complicated and interesting in contrast to 
the simple structure of the boundary layer for slow 
compaction. Since A > 1 and the highly nonlin- 
ear permeability function k = (</>/0 o ) m , m = 8, 
the governing equations are also highly nonlinear. 
However, we use these features and pursue asymp- 
totic analysis to seek appropriate asymptotic solu- 
tions. 

4.3 Poroelastic Compaction 

For the case A ^> 1, we rewrite (25) as 

£+«£n-| -<!-*>]£ --(.-«, 
-(T^(!+<n-I-< 1 -«iI>.« 

By using the perturbations 
+ 



A 



-rhW + ..., P = p(°) + ip « + ..., (41) 
A 



the leading order equation becomes 
90(0) ^(o) Qp(o) 



(42) 



dz (1-0W) dz ' 

whose integration gives 

p(°)=ln(0 o /0( o ))-(0o-0 (O) ), (43) 

which is the Athy-type relation and is exactly the 
same form as in Smith (1971) and Fowler and Yang 
(1998). From (24) the equation of leading order is 
thus 



or 



(l-0(°))d0(°) 
0(0) 

g0_W 
dz ~ 



(1-0(°)) = O, (44) 



0. 



The boundary condition 0(°) = O gives 



0(°) 



b e 



-(h-z) 



(45) 



(46) 



which decreases with depth h — z exponentially. 
This solution is the same as the equilibrium solu- 
tion in the poroelastic case, and thus the top region 
of viscoelastic compaction is essentially poroelastic 
and the viscous effect is only of secondary impor- 
tance in this region. However, as (f> decreases, the 
term A(0/</>o) m may become very small due to the 



large exponent m = 8. The relation A(0/0o) m = 1 
defines a critical value of <j> in the transition region 



i) e 



(47) 



In fact, the above solutions are only valid when 
0(°) > <p* and h - z < n = (In A)/m. 

4.4 Transition Region 

Since <fi <~ 0* , we define 



= 0*e 



^-n + ^i^, (48) 



m 



— > P^P , 

to to 



(49) 



wherep* = ln(0o/0*) — (<t>o— 0*)- By changing vari- 
ables to (t, T]) via dt — > 9t — mhd n , d z — > m9^, and 
assuming to ^ 1, we have the equation of leading 
order 

+ (1 - 0*)W„ = 0, (50) 
W = e*[P v -{l-<F)], (51) 



Thus 



(1-0*)" 



W = W* 



hcj)* 



(1-0*) 



(53) 



^ - 1 = (1 f )p * + (^ - ^) e -* (54) 



ft,0* 



-0o 



(1-0*)!^* 

hcj)* 



(55) 



whose solution can be written as a quadrature. 
In the limit 77 — > —00, P v — > 0, equation (52) 
shows that the dominant term is the viscous term 
(1 — 4>*)p* so that compaction will gradually shift 
from viscoelastic to purely viscous behavior. This 
has important geophysical implication for com- 
paction in sedimentary basins, since the purely vis- 
cous mechanism may be responsible for overpres- 
suring and mineralized seals in oil-reservoir and 
hydrocarbon basins. Furthermore, In order to de- 
termine h, we require (53) and (55) to match the 
solution below the transition layer. 

4.5 Viscous Compaction 

In the region below the transition layer, the poros- 
ity < 0* is usually very small, while the effective 
pressure p is increasing andp <~ p* = 0(1). Rewrit- 
ing (25) as 



du s 
dz 



dp 
[ 8t 



s dp } 



(56) 



From (49) and (52), we know that p changes slowly, 
which implies § ~ §E « 1 or 0(§f + U *|E ) « p . 
We then have approximately 



p w -■ 



du s 
dz 



(57) 



which implies that compaction is now essentially 
purely viscous. Thus we get 



! = (58) 
6 d 2 u s 



u = 



which are the equations solved by Fowler and Yang 
(1999) when S = 1 for purely viscous compaction. 
Following the same solution procedure given by 
Fowler and Yang (1999), we can expect to get the 
same solutions. Thus, we only write down here the 
solution for h 



L ■ r(Wo) 

ft = TO, 1 



(1-0*) to(1 



2p* 



In to, (60) 



where 7 = £ ■ This essentially completes 

the solution procedure. Figure 2 shows the compar- 
ison of numerical results with the above obtained 
asymptotic solutions (46) and (55) in the poroelas- 
tic and transition region. 

5 Discussion 

The present model of viscoelastic flow and nonlin- 
ear compaction in sedimentary basins uses a rheo- 
logical relation which incorporates both poroelastic 
and viscous effects in the 1-D compacting frame. 
Based on the frame invariance of irrotational fea- 
ture of the 1-D flow, a generalised viscoelastic com- 
paction relation of Maxwell type has been formu- 
lated. The nondimensional model equations are 
mainly controlled by one parameter A, which is the 
ratio of the hydraulic conductivity to the sedimen- 
tation rate. Following a similar asymptotic analysis 
given by Fowler and Yang (1998) , we have been able 
to obtain the approximate solutions for either slow 
compaction (A <C 1) or fast compaction (A ^> 1). 
The more realistic and yet more interesting case is 
when ) » 1, and the solution gives a nearly ex- 
ponential profile of porosity versus depth, which 
implies that compaction in the top region is es- 
sentially poroelastic and its profile is virtually at 
equilibrium. 

The numerical simulations and asymptotic anal- 
ysis have shown that porosity-depth profile is 




on pressure solution and its application to some 
field problems such as land subsidence associ- 
ated with fluid withdrawal from undercompacted 
aquifers, Revil (1999) suggests a Voigt-type poro- 
visco-plastic rheological behavior to characterize 
pressure solution and to applications to some field 
problems including equilibrium and disequilibrium 
compactions and subsidence. Naturally, more work 
- is needed to incorporate a Voigt-type rheology ap- 
plied to compaction in addition to the present 
Maxwell-type law. 
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Figure 2: Numerical results (solid curves) with 
A = 100, t = 10. The asymptotic solutions (46) 
and (55) are also plotted as a comparison (dashed). 
Profile in the top region is nearly exponential fol- 
lows by a transition to pure viscous compaction 
where porosity is nearly uniform. 
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